This is an R Markdown document containing code for the workshop “Accessing the USGS National Map and making 3D maps with terrainr”, held virtually on 2021-05-28. If you want to follow along with the workshop as we go, click this link to download the R Markdown notebook.

If you’re not familiar with R Markdown, this document lets us write both plain text and code in a single document. Document sections inside three ` marks are called “code chunks”:

plot(Orange)

A single line in a code chunk can be run by pressing Control+Enter with your cursor on the line. The entire chunk can be run by pressing Control+Shift+Enter with your cursor inside the chunk.

Part 1: Installing Packages

For this workshop, we first need to install the packages that we’ll be using. This chunk will install any missing packages, and then load them all:

if (!require("sf")) install.packages("sf")
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
if (!require("terrainr")) install.packages("terrainr")
## Loading required package: terrainr
if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
if (!require("raster")) install.packages("raster")
## Loading required package: raster
## Loading required package: sp
if (!require("progressr")) install.packages("progressr")
## Loading required package: progressr
if (!require("progress")) install.packages("progress")
## Loading required package: progress
library("sf")
library("terrainr")
library("ggplot2")
library("raster")
library("progressr")
handlers("progress")

Part 2: Retrieving Data

Up next, let’s download some data! For the purposes of today’s workshop, we’ll work with campsite locations within Bryce Canyon National Park. We can download this data from the National Park Service using the read_sf function from sf:

campsites <- read_sf("https://opendata.arcgis.com/datasets/06993a2f20bf42d382c8ce6bd54027c4_0.geojson")

This dataset can be found online at https://public-nps.opendata.arcgis.com/datasets/bryce-canyon-national-park-backcountry-campsites-?geometry=-112.494%2C37.471%2C-111.908%2C37.566

We now have an sf object containing a handful of campsite locations! We can plot our data using ggplot2 to get a sense of what the spatial distribution of the campsites looks like:

ggplot(campsites) + 
  geom_sf(shape = 4, color = "red") + 
  theme_void()

This is the area we want to download data for from the National Map. Note that everything we do from here on out works with pretty much any sf object; for instance, the terrainr documentation (at https://docs.ropensci.org/terrainr/) often starts with a single latitiude/longitude point and then works from there.

By default, terrainr will download data for the bounding box of whatever sf object you provide it. That is, data will only be downloaded for the smallest rectangle of area that includes your entire sf object.

In this case, we want a little bit of a buffer around our campsite locations so that our points aren’t at the exact edge of the map. To make that happen, we can use set_bbox_side_length to specifically control how much data to download. In this case, a 16 kilometer side length works pretty well:

campsite_bbox <- set_bbox_side_length(campsites, 16, "km")

We can add this to our plot to make sure we’re giving a good buffer to our furthest points:

ggplot() + 
  geom_sf(data = st_set_crs(campsite_bbox, 4326), color = "black", fill = NA) + 
  geom_sf(data = campsites, shape = 4, color = "red") + 
  theme_void()

Since all of our points fall comfortably within this box, we can use it to download our data. We’ll use the get_tiles function to query the National Map to download elevation data (from the USGS 3D Elevation Program DEM) and orthoimagery (from the USDA National Agricultural Imagery Program). This function will break our query into smaller parts and then save the tiles as separate files; we can control where the files get saved by using the output_prefix argument.

We’ll wrap this query in the with_progress function in order to have a progress bar display while our data downloads. While for the purposes of this workshop we’ll be downloading 30 meter rasters, which download pretty quickly (and only require one tile per type of data), this can be very useful for monitoring downloads when using higher-resolution data.

with_progress(
  output_tiles <- get_tiles(campsite_bbox, 
                            output_prefix = "bryce_canyon_np_30m",
                            services = c("elevation", "ortho"),
                            resolution = 30)
)

By assigning the output from get_tiles to output_tiles, we create a list of all of the map tiles we’ve downloaded:

output_tiles
## $elevation
## [1] "bryce_canyon_np_30m_3DEPElevation_1_1.tif"
## 
## $ortho
## [1] "bryce_canyon_np_30m_USGSNAIPPlus_1_1.tif"

Part 3: Plotting Data

This can be really helpful when trying to work with these files inside R. For instance, to see our elevation raster, we can load it into R using the raster function, then graph it using plot:

elevation <- raster(output_tiles$elevation)

plot(elevation, col = terrain.colors(255))

We can visualize our orthoimagery in a similar way. Because this image is stored as a multi-band raster (that is, as a color image), we’ll have to use the stack and plotRGB functions in the place of raster and plot:

ortho <- stack(output_tiles$ortho)

plotRGB(ortho, scale = 1)

These visualizations are fast and work perfectly fine, but can be difficult to extend further. For a bit more flexibility, we might try and use the ggplot2 package to plot our data instead. Plotting the elevation raster can be done easily enough using functions from ggplot2. First, we’ll convert the raster to a data frame:

elevation <- as.data.frame(elevation, xy = TRUE)
names(elevation) <- c("x", "y", "z")

And then plot the data using geom_raster:

ggplot(elevation, aes(x, y, fill = z)) +
  geom_raster() + 
  scale_fill_gradientn(colors = terrain.colors(255)) + 
  coord_sf(crs = 4326) + 
  theme_void()

Plotting our orthoimagery is again a bit of a different process. This time, we need to use a new function from terrainr called geom_spatial_rgb:

ggplot() + 
  geom_spatial_rgb(data = ortho,
                   mapping = aes(x, y, r = red, g = green, b = blue)) + 
  coord_sf(crs = 4326) + 
  theme_void()

Notice we don’t need to convert our data to a data frame, like we did with elevation. When plotting basemaps with geom_spatial_rgb, you can provide RasterStack objects created by stack, or data frames, or file paths to the images you want to load.

Using ggplot2 makes it a lot easier to, for instance, plot our campsites on top of a base map:

ggplot() + 
  geom_spatial_rgb(data = ortho,
                   mapping = aes(x, y, r = red, g = green, b = blue)) + 
  geom_sf(data = campsites, shape = 4, color = "red") + 
  coord_sf(crs = 4326) + 
  theme_void()

We can then continue to edit this map just like we would any other ggplot.

terrainr also provides a few options for editing our basemaps. If we want to, for instance, add our campsites to the orthoimagery, we can first create a new image from our campsite locations:

vector_to_overlay(
  campsites,
  output_tiles$ortho,
  "campsite_overlay.png",
  shape = 4,
  color = "red"
)

This new basemap can be loaded and displayed just like our orthoimage map:

campsite_overlay <- stack("campsite_overlay.png")
plotRGB(campsite_overlay)
## Warning in .local(x, ...): layer was changed to 1

## Warning in .local(x, ...): layer was changed to 1

To combine this new image with our orthoimagery, we’ll want to use the combine_overlays function:

combine_overlays(
  output_tiles$ortho,
  "campsite_overlay.png",
  output_file = "campsites_ortho.png"
)

This function “stacks” images, putting the first image in the set at the bottom of the pile and the last image on the top. As a result, we now have our campsite locations on top of our orthoimagery, as a single image:

campsites_ortho <- stack("campsites_ortho.png")
plotRGB(campsites_ortho)

This process is particularly useful for combining multiple map tiles from the National Map into a single basemap, allowing you to overlay roads, rivers, and more on orthoimagery.

Part 4: Unity

This part of the workshop will be performed as a demonstration, rather than a tutorial, due to the amount of time it takes to download the data while on video calls. The code used in the Unity demonstration is below – fair warning, it takes a (very) long time to run!

with_progress(
  output_tiles <- get_tiles(campsite_bbox, 
                            output_prefix = "bryce_canyon_np",
                            services = c("elevation", "ortho"))
)

mapply(
  merge_rasters,
  input_rasters = output_tiles,
  output_raster = paste0(names(output_tiles), ".tif")
)

mapply(
  raster_to_raw_tiles,
  input_file = c("elevation.tif", "ortho.tif"),
  output_prefix = "bryce_canyon",
  raw = c(TRUE, FALSE)
)